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Abstract 

We present a mathematical analysis of a networks with Integrate-and-Fire neurons with con- 
ductance based synapses. Taking into account the realistic fact that the spike time is only known 
within some finite precision, we propose a model where spikes are effective at times multiple of a 
characteristic time scale S, where S can be arbitrary small (in particular, well beyond the numerical 
precision). We make a complete mathematical characterization of the model-dynamics and obtain 
the following results. The asymptotic dynamics is composed by finitely many stable periodic orbits, 
whose number and period can be arbitrary large and can diverge in a region of the synaptic weights 
space, traditionally called the "edge of chaos", a notion mathematically well defined in the present 
paper. Furthermore, except at the edge of chaos, there is a one-to-one correspondence between the 
membrane potential trajectories and the raster plot. This shows that the neural code is entirely "in 
the spikes" in this case. As a key tool, we introduce an order parameter, easy to compute numeri- 
cally, and closely related to a natural notion of entropy, providing a relevant characterization of the 
computational capabilities of the network. This allows us to compare the computational capabilities 
of leaky and Integrate-and-Fire models and conductance based models. The present study considers 
networks with constant input, and without time-dependent plasticity, but the framework has been 
designed for both extensions. 
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Neuronal networks have the capacity to treat incoming information, performing complex computa- 
tional tasks (see (Rieke, Warland, Steveninck, & Bialek, 1996) for a deep review), including sensory-motor 
tasks. It is a crucial challenge to understand how this information is encoded and transformed. However, 
when considering in vivo neuronal networks, information treatment proceeds usually from the interaction 
of many different functional units having different structures and roles, and interacting in a complex 
way. As a result, many time and space scales are involved. Also, in vivo neuronal systems are not 
isolated objects and have strong interactions with the external world, that hinder the study of a specific 
mechanism (Fregnac, 2004). In vitro preparations are less subject to these restrictions, but it is still 
difficult to design specific neuronal structure in order to investigate the role of such systems regarding 
information treatment (Koch & Segev, 1998). In this context models are often proposed, sufficiently close 
from neuronal networks to keep essential biological features, but also sufficiently simplified to achieve a 
characterization of their dynamics, the most often numerically and, when possible, analytically (Gerstner 
& Kistler, 2002b; Dayan & Abbott, 2001). This is always a delicate compromise. At one extreme, one 
reproduces all known features of ionic channels, neurons, synapses ... and lose the hope to have any 
(mathematics and even numeric) control on what is going on. At the other extreme, over-simplified mod- 
els can lose important biological features. Moreover, sharp simplifications may reveal exotic properties 
which are in fact induced by the model itself, but do not exist in the real system. This is a crucial aspect 
in theoretical neuroscience, where one must not forget that models are subject to hypothesis and have 
therefore intrinsic limits. 

For example, it is widely believed that one of the major advantages of the Integrate-and-Fire (IF) 
model is its conceptual simplicity and analytical tractability that can be used to explore some general 
principles of neurodynamics and coding. However, though the first IF model was introduced in 1907 
by Lapicquc (Lapicque, 1907) and though many important analytical and rigorous results have been 
published, there are essential parts missing in the state of the art in theory concerning the dynamics 
of IF neurons (see e.g. (Mirollo & Strogatz, 1990; Ernst, Pawelzik, & Gcisel, 1995; Senn & Urbanczik, 
2001; Timme, Wolf, & Geisel, 2002; Memmeshcimer & Timmc, 2006; Gong & Leeuwcn, 2007; Jahnke, 
Memmesheimer, & Timme, 2008) and references below for analytically solvable network models of spiking 
neurons). Moreover, while the analysis of an isolated neuron submitted to constant inputs is straight- 
forward, the action of a periodic current on a neuron reveals already an astonishing complexity and the 
mathematical analysis requires elaborated methods from dynamical systems theory (Keener, Hoppen- 
steadt, & Rinzel, 1981; Coombes, 1999b; Coombes & Bressloff, 1999). In the same way, the computation 
of the spike train probability distribution resulting from the action of a Brownian noise on an IF neuron is 
not a completely straightforward exercise (Knight, 1972; Gerstner & Kistler, 2002a; Brunei & Sergi, 1998; 
Brunei & Latham, 2003; Touboul & Faugeras, 2007) and may require rather elaborated mathematics. At 
the level of networks the situation is even worse, and the techniques used for the analysis of a single neuron 
are not easily extensible to the network case. For example, Bressloff and Coombes (Bressloff & Coombes, 
2000b) have extended the analysis in (Keener et al., 1981; Coombes, 1999b; Coombes & Bressloff, 1999) 
to the dynamics of strongly coupled spiking neurons, but restricted to networks with specific architectures 
and under restrictive assumptions on the firing times. Chow and Kopell (Chow & Kopell, 2000) studied 
IF neurons coupled with gap junctions but the analysis for large networks assumes constant synaptic 
weights. Brunei and Hakim (Brunei & Hakim, 1999) extended the Fokkcr-Planck analysis combined to 
a mean-field approach to the case of a network with inhibitory synaptic couplings but under the as- 
sumptions that all synaptic weights are equal. However, synaptic weight variability plays a crucial role 
in the dynamics, as revealed e.g. using mean-field methods or numerical simulations (VanVrceswijk & 
Hansel, 1997; VanVreeswijk & Sompolinsky, 1998; VanVreeswijk, 2004). Mean-field methods allow the 
analysis of networks with random synaptic weights (Amari, 1972; Sompolinsky, Crisanti, & Sommers, 
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1988; Cessac, Doyon, Quoy, & Samuclidcs, 1994; Cessac, 1995; Hansel & Mato, 2003; Soula, Beslon, & 
Mazet, 2006; Samuelides & Cessac, 2007) but they require a "thermodynamic limit" where the number of 
neurons tends to infinity and finite-size corrections are rather difficult to obtain. Moreover, the rigorous 
derivation of the mean-field equations, that requires large-deviations techniques (BenArous & Guionnet, 
1995), has not been yet done for the case of IF networks with continuous time dynamics (for the discrete 
time case see (Soula et al., 2006; Samuelides & Cessac, 2007)). 

Therefore, the "analytical tractability" of IF models is far from being evident. In the same way, the 
"conceptual simplicity" hides real difficulties which are mainly due to the following reasons. IF models 
introduce a discontinuity in the dynamics whenever a neuron crosses a threshold: this discontinuity, 
that mimics a " spike" , maps instantaneously the membrane potential from the threshold value to a reset 
value. The conjunction of continuous time dynamics and instantaneous reset leads to real conceptual 
and mathematical difficulties. For example, an IF neuron without refractory period (many authors have 
considered this situation), can, depending on parameters such as synaptic weights, fire uncountably 
many spikes within a finite time interval, leading to events which are not measurable (in the sense of 
probability theory). This prevents the use of standard methods in probability theory and notations 
such as p(t) — Y^ii=i ^ — ^*) (spike response function ) simply lose their meaning 1 . Note also that the 
information theory (e.g. the Shannon theorem, stating that the sampling period must be less than half 
the period corresponding to the highest signal frequency) is not applicable with unbounded frequencies. 
But IF models have an unbounded frequencies spectrum (corresponding to instantaneous reset). From 
the information theoretic point of view, it is a temptation to relate this spurious property to the erroneous 
fact that the neuronal network information is not bounded. These few examples illustrate that one should 
not be abused by the apparent simplicity of IF models and must be careful in pushing too much their 
validity in order to explore some general principles of neurodynamics and coding. 

The situation is not necessarily better when considering numerical implementations of IF neurons. 
Indeed, it is known from a long time that the way the membrane potential is reset in a neuronal network 
simulation have significant consequences for the dynamics of the model. In particular, Hansel et al 
(Hansel, Mato, Meunier, & Neltner, 1998) showed that a naive implementation of integrate-and-fire 
dynamics on a discrete time grid introduces spurious effects and proposed an heuristic method to reduce 
the errors induced by time discretization. In parallel, many people have developed event based integration 
schemes (Brette et al., 2007), using the fact that the time of spike of a neuron receiving instantaneous 
spikes from other neurons can be computed analytically, thus reducing consequently the computation 
time and affording the simulation of very large networks. In addition, exact event based computational 
schemes are typically used for the above-mentioned analytically tractable model classes, see, e.g. (Mirollo 
& Strogatz, 1990; Timme et al., 2002). Unfortunately, this approach suffers two handicaps. If one 
considers more elaborated models than analytically tractable models, one is rapidly faced to the difficulty 
of finding an analytical expression for the next spike time (Rudolph & Dcstexhe, 2006). Moreover, any 
numerical implementations of a neural network model will necessarily introduce errors compared to the 
exact solution. The question is : how does this error behave when iterating the dynamics ? Is it amplified 
or damped ? In IF models, as set previously, these errors are due to the discontinuity in the membrane 
potential reset and to the time discretization. This has been nicely discussed by Hansel et al (Hansel et 
al., 1998) paper. These authors point out two important effects. When a neuron fires a spike between 
time t and t + At a local error on the firing time is made when using time discretization. First, this leads 
to an error on the membrane potential and second this error is propagated to the other neurons via the 

1 Obviously, one can immediately point out that (i) this situation is not plausible if one thinks of biological neurons 
and (ii) is not "generic" for IF models. Thus, objection (i) implies that some conclusions drawn from IF models are not 
biologically plausible, while objection (ii) needs to be made mathematically clear. This is one of the goals of this paper. 
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synaptic interaction term. Unfortunately, this analysis, based on numerical simulations, was restricted to 
a specific architecture (identical excitatory neurons) and therefore the conclusions drawn by the authors 
cannot be extended as it is to arbitrary neural architectures. Indeed, as we show in the present paper, 
the small error induced by time discretization can be amplified or damped, depending on the synaptic 
weights value. This leads to the necessity of considering carefully (that is mathematically) the spurious 
effects induced by continuous time and instantaneous reset in IF models, as well as the effects of time 
discretization. This is one aspect discussed in the present paper. 

More generally, this work contains several conclusions forming a logical chain. After a discussion on 
the characteristic times involved in real neurons and comparison to the assumptions used in Integrate and 
Fire models we argue that discrete time IF models with synchronous dynamics can be used to model real 
neurons as well, provided that the time scale discretization is sufficiently small. More precisely, we claim 
that IF equations are inappropriate if one sticks to much on the instantaneous reset and spike time, but 
that they provide a good and mathematically tractable model if one allows reset and spike to have some 
duration. We therefore modify the reset and spike definition (while keeping the differential equation for 
the dynamics of the membrane potential below the threshold). The goal is however NOT to propose yet 
another numerical scheme for the numerical integration of continuous time IF models. Instead, our aim 
is to analyze mathematically the main properties of the corresponding dynamical system, describing the 
evolution of a network with an arbitrary, finite, size (i.e. we don't use neither a mean-field approach nor 
a thermodynamic limit). We also consider an arbitrary architecture. Finally, in our analysis the time 
discretization step is arbitrary small, (thus possibly well below the numerical precision). For this, we 
use a dynamical system approach developed formerly in (Blanchard, Cessac, & Krueger, 2000; Cessac, 
Blanchard, Krueger, & Meunier, 2004). In particular, in (Cessac, 2008), a discrete time version of a leaky 
Integrate-and-Fire network, was studied. It was shown that the dynamics is generically periodic, but the 
periods can become arbitrary large (in particular, they can be larger than any accessible computational 
time) and in (non generic) regions of the synaptic weights space, dynamics is chaotic. In fact, a complete 
classification of the dynamical regimes exhibited by this class of IF models was proposed and a one-to- 
one correspondence between membrane potential trajectories and raster plots was exhibited (for recent 
contributions that study periodic orbits in large networks of IF neurons, see (Jahnkc et al., 2008; Gong & 
Leeuwen, 2007).). Beyond these mathematical results, this work warns one about some conclusions drawn 
from numerical simulations and emphasizes the necessity to have, when possible, a rigorous analysis of 
the dynamics. 

The paper (Cessac, 2008) dealt however with a rather simple version of IF neurons (leaky Integrate 
and Fire) and one may wonder whether this analysis extend to models closer to biology. In the present 
paper we extend these results, and give a mathematical treatment of the dynamics of spikes generated 
in synaptic coupled integrate-and-fire networks where synaptic currents are modeled in a biophysically 
plausible way (conductance based synapses). As developed in the text, this extension is far from being 
straightforward and requires a careful definition of dynamics incorporating the integration on the spikes 
arising in the past. This requires a relatively technical construction but this provides a setting where a 
rigorous classification of dynamics arising in IF neural networks with conductance based synapse can be 
made, with possible further extended to more elaborated models. 

The paper is organized as follows. In section 1 we give a short presentation of continuous time 
Integrate-and-Fire models. Then, a careful discussion about the natural time scales involved in biological 
neurons dynamics and how continuous time IF models violate these conditions is presented. From this 
discussion we propose the related discrete time model. Section 2 makes the mathematical analysis of 
the model and mathematical results characterizing its dynamics are presented. Moreover, we introduce 
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an order parameter, called d(f2,5), which measures how close to the threshold are neurons during their 
evolution. Dynamics is periodic whenever d(Ct, S) is positive, but the typical orbit period can diverge 
when it tends to zero. This parameter is therefore related to an effective entropy within a finite time 
horizon, and to the neural network capability of producing distinct spikes trains. In other words, this is 
a way to measure the ability of the system to emulate different input-output functions. See (Langton, 
1990; Bertschinger & Natschlager, 2004) for a discussion on the link between the system dynamics 
and its related computational complexity 2 . The smaller d(f2,<S), the larger is the set of distinct spikes 
trains that the neural network is able to produce. This implies in particular a larger variability in the 
responses to stimuli. The vanishing of d(Q, S) corresponds to a region in the parameters space, called 
"the edge of chaos" , and defined here in mathematically precise way. In section 3 we perform numerical 
investigations of d($l,S) in different models from leaky Integrate-and-Fire to conductance based models. 
These simulations suggest that there is a wide region of synaptic conductances where conductance based 
models display a large effective entropy, while this region is thinner for leaky Integrate-and-Fire models. 
This provides a quantitative way to measuring how conductances based synapses and currents enhances 
the information capacity of integrate and fire models. Section 4 proposes a discussion on these results. 

1 General framework. 

1.1 General structure of Integrate and Fire models. 

We consider the (deterministic) evolution of a set of N neurons. Call Vk(t) the membrane potential of 
neuron k € {1 . . . N} at time t and let V(£) be the vector [Vk(t)]^ =1 . We denote by V = V(0) the initial 
condition and the (forward) trajectory of V by: 

v {v(t)}£ , 

where time can be either continuous or discrete. In the examples considered here the membrane potential 
of all neurons is uniformly bounded, from above and below, by some values V m in,V m ax- Call M = 
[Vmin, V m ax] N ■ This is the phase space of our dynamical system. 

We are focusing here on "Integrate and Fire models" , which always incorporate two regimes. For the 
clarity of the subsequent developments we briefly review these regimes (in a reverse order) . 

1. The "fire" regime. 

Fix a real number 9 6 \V m i n ,V ma x] called the firing threshold of the neuron 3 . Define the firing 

2 It has been proposed that optimal computational capabilities are achieved by systems whose dynamics is neither 
chaotic nor ordered but somewhere in-between order and chaos. This has led to the idea of computation at "the edge 
of chaos". Early evidence for this hypothesis has been reported by Kauffman (Kauffman, 1969) and Langton (Langton, 
1990) considering cellular automata behavior, and Packard (Packard, 1988) using a genetic algorithm. See (Bertschinger 
& Natschlager, 2004) for a review. In relation, with these works, theoretical results by Derrida and co-authors (Dcrrida 
& Pomeau, 1986; Derrida & Flyvbjerg, 1986) allow to characterize analytically the dynamics of random Boolean networks 
and for networks of threshold elements (Derrida, 1987). Recently (Bertschinger & Natschlager, 2004) have contributed to 
this question, considering numerical experiments in the context of real-time computation with recurrent neural networks. 

3 We assume that all neurons have the same firing threshold. The notion of threshold is already an approximation 
which is not sharply defined in Hodgkin-Huxley (Hodgkin & Huxley, 1952) or Fitzhugh-Nagumo (FitzHugh, 1961; Nagumo, 
Arimoto, & Yoshizawa, 1962) models (more precisely it is not a constant but it depends on the dynamical variables). Recent 
experiments (Naundorf, Wolf, & Volgushev, 2006; McCormick, Shu, & Yu, 2007; Naundorf, Wolf, & Volgushev, 2007) even 
suggest that there may be no real potential threshold. 
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times of neuron fc, for the trajectory 4 V, by: 

4 n) (V) = mf{t \t > tjrV), V k {t) > 9} (1) 

where t k °^ = —00. The firing of neuron k corresponds to the following procedure. If Vk(t) > 9 then 
neuron membrane potential is reset instantaneously to some constant reset value V reset and a spike 
is emitted toward post-synaptic neurons. In mathematical terms firing reads 5 : 

V k (t) >9^V k (t+) = V reset (2) 

where V reset <E [V m i n , V max ] is called the "reset potential". In the sequel we assume, without loss 
of generality, that V rese t — 0. This reset has a dramatic effect. Changing the initial values of 
the membrane potential, one may expect some variability in the evolution. Now, fix a neuron k 
and assume that there is a time t > and an interval [a, b] such that, VVfe(O) G [a, b], V k (t) > 9. 
Then, after reset, this interval is mapped to the point V reset - Then, all trajectories born from [a, b] 
collapse on the same point and have obviously the same further evolution. Moreover, after reset, 
the membrane potential evolution does not depend on its past value. This induces an interesting 
property used in all the Integrate and Fire models that we know (see e.g. (Gerstner & Kistler, 
2002b)). The dynamical evolution is essentially determined by the firing times of the neurons, 
instead of their membrane potential value. 

2. The "Integrate regime". 

Below the threshold, V k < 9, neuron fc's dynamics is driven by an equation of form: 

C^+g k V k = i k , (3) 

where C is the membrane capacity of neuron k. Without loss of generality we normalize the 
quantities and fix C = 1. In its most general form, the neuron fc's membrane conductance g k > 
depends on V k (see e.g. Hodgkin-Huxley equations (Hodgkin & Huxley, 1952)) and time t, while 
the current i k can also depend on V, the membrane potential vector, on time t, and also on the 
collection of past firing times. The current i k can include various phenomenological terms. Note 
that (3) deals with neurons considered as points instead of spatially extended objects. 

Let us give two examples investigated in this paper. 

(i) The leaky integrate and fire model. 

In its simplest form equation (3) reads: 

4 Note that, since the dynamics is deterministic, it is equivalent to fix the forward trajectory or the initial condition 
V = V(0). 

5 Note that the firing condition includes the possibility to have a membrane potential value above the threshold. This 
extension of the standard definition affords some discontinuous jumps in the dynamics. These jumps arise when considering 
addition of (discontinuous) noise, or a profiles with jumps (e.g. a(t) = i-e~ ~ , t > 0). They also appear when considering 
a discrete time evolution. Note that strictly speaking, this can happen, within the numerical precision, even with numerical 
schemes using interpolations to locate more precisely the spike time (Hansel et al., 1998). 
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(4) 



where g k is a constant, and r k = |j is the characteristic time for membrane potential decay when 
no current is present. This model has been introduced in (Lapicque, 1907). 

(ii) Conductance based models with a profiles 

More generally, conductance and currents depend on V only via the previous firing times of the 
neurons (Rudolph & Destexhe, 2006). Namely, conductances (and currents) have the general form 6 , 

9k = 9k {t, I ) where is the n-th firing time of neuron j and | is the list of firing times 
of all neurons up to time t. This corresponds to the fact that the occurrence of a post-synaptic 
potential on synapse j, at time t^ n \ results in a change of the conductance g k of neuron k. As an 
example, we consider models of form: 

^ = ~ (Vk E L ) - i syn \ Vk ,t, {tf} t ) + *<»«(*) (5) 
where the first term in the r.h.s. is a leak term, and where the synaptic current reads: 

N N 

i^ n \V k ,t, {tf}) = (V k E+) J>+ (t, {tf}) + (V k E-) £^(t, {tf} t ), 
where E ± are reversal potential (typically E + ~ OmV and E~ ~ — 75mV) and where: 

M^t.V) 

In this equation, Mj(t, V) is the number 7 of times neuron j has fired at time t. G^- is the synaptic 
efficiency (or synaptic weight) of the synapse j — ► k. (It is zero if there is no synapse j — > k), 
where + [— ] expresses that synapse j — > k is excitatory [inhibitory]. The a function mimics the 
conductance time-course after the arrival of a post-synaptic potential. A possible choice is: 

a ± (t) = H(t)^e-^, (6) 



6 Thc rather cumbersome notation g k {t, j^'j™' ) simply expresses that in conductance based models the conductance 

depends on the whole set (history) of (past) firing times. Note that membrane potentials are reset after neuron firing, but 
not neuron conductances. 

^Henceforth, one assumes that there are finitely many spikes within a finite time interval. For continuous time dynamics, 

this fact is not guaranteed when neglecting the refractory period. Note also that this number, as well as the list j^j™' , 

depends on the initial condition V and a small change in the initial condition may induce a drastic change of Mj (t, V) at 
time t, as discussed later. This effect is sometimes disregarded (Coombcs, 1999b). This issue has also been discussed (for 
current based IF-likc models) as "phase history functions" in (Ashwin & Timmc, 2005; Broer, Efstathiou, & Subramanian, 
2008) (we thank one of the reviewers for this remark). 
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with H the Heaviside function and r being characteristic times. This synaptic profile, with 
a(0) = while a(t) is maximal for t = r, allows us to smoothly delay the spike action on the 
post-synaptic neuron. We are going to neglect other forms of delays in the sequel. 

Then, we may write (5) in the form (3) with: 

N Mj(t,V) N M 3 (t,V) 

*(*.{*S B) } t ) = ^+£<% E « + M n) )+E^- E «-(*-«<->), (7) 

j — 1 n— 1 j — 1 n— 1 

and: 

N N 

i k (t, {t^}) = ^+E+ 5>+(*. {*W} t ) + ZT £^( t) {^} t ) + £(«*)(«). (8) 
1.2 Discrete time dynamics. 

1.2.1 Characteristic time scales in neurons dynamics. 

Integrate and Fire models assume an instantaneous reset of the membrane potential corollary to an 
infinite precision for the spike time. We would like to discuss shortly this aspect. Looking at the spike 
shape reveals some natural time scales: the spike duration r (a few ms); the refractory period r ~ 1ms; 
and the spike time precision. Indeed, one can mathematically define the spike time as the time where the 
action potential reaches some value (a threshold, or the maximum of the membrane potential during the 
spike), but, on practical ground, spike time is not determined with an infinite precision. An immediate 
conclusion is that it is not correct, from an operational point of view, to speak about the "spike time" , 
unless one precises that this time is known with a finite precision St. Thus the notion of list of firing 

time 1 4"^} us °d m section 1.1, must be revisited, and a related question is "what is the effect of this 
indeterminacy on the dynamical evolution ?" . Note that this (evident ?) fact is forgotten when modeling 
e.g. spike with Dirac distributions. This is harmless as soon as the characteristic time St is smaller 
than all other characteristic times involved in the neural network. This is essentially true in biological 
networks but it is not true in Integrate and Fire models. 

These time scales arise when considering experimental data on spikes. When dealing with models, 
where membrane potential dynamics is represented by ordinary differential equations usually derived 
from Hodgkin-Huxley model, other implicit times scales must be considered. Indeed, Hodgkin-Huxley 
formulation in term of ionic channel activity assumes an integration over a time scale dt which has to be 
(i) quite larger than the characteristic time scale Tp of opening/closing of the channels, ensuring that the 
notion of probability as a meaning; (ii) quite larger than the correlation time tq between channel states 
ensuring that the Markov approximation used in the equations of the variable to, n, h is legal. This means 
that, although the mathematical definition of assumes a limit dt — > 0, there is a time scale below which 
the ordinary differential equations lose their meaning. Actually, the mere notion of "membrane potential" 
already assumes an average over microscopic time and space scales. Note that the same is true for all 
differential equations in physics ! But this (evident ?) fact is sometimes forgotten when dealing with 
Integrate and Fire models. Indeed, to summarize, the range of validity of an ODE modeling membrane 
potential dynamics is max(Tc,Tp) << dt « St << r. But the notion of instantaneous reset implies 
t = and the mere notion of spike time implies that St — !! 
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There is a last time scale related to the notion of raster plot. It is widely admitted that the "neural 
code" is contained in the spike trains. Spike trains are represented by raster plots, namely bi-dimensional 
diagrams with time on abscissa and some neurons labeling on ordinate. If neuron k fires a spike "at 
time ifc" one represents a vertical bar at the point (tk,k). Beyond the discussion above on the spike 
time precision, the physical measurement of a raster plot involves a time discretization corresponding to 
the time resolution 5a of the apparatus. When observing a set of neurons activity, this introduces an 
apparent synchronization, since neurons firing between t and t + 5 a will be considered as firing simul- 
taneously. This raises several deep questions. In such circumstances the "information" contained in the 
observed raster plot depends on the time resolution 5a (Golomb., Hertz, Panzeri, Treves, & Richmond, 
1997; Panzeri & Treves, 1996) and it should increase as 5a decreases. But is there a limit time resolution 
below which this information does not grow anymore ? In Integrate and Fire models this limit is 5a = 0. 
This may lead to the conclusion that neural networks have an unbounded information capacity. But is 
this a property of real neurons or only of Integrate and Fire models ? 

The observation of raster plots corresponds to switching from the continuous time dynamics of mem- 
brane potential to the discrete time and synchronous dynamics of spike trains. One obtains then, in 
some sense, a new dynamical system, of symbolic type, where variables are bits ('0' for no spike, and T' 
otherwise) . The main advantage of this new dynamical system is that it focuses on the relevant variables 
as far as information and neural coding is concerned i.e. one focuses on spikes dynamics instead of mem- 
brane potentials. In particular, membrane potentials may still depend continuously on time, but one is 
only interested in their values at the times corresponding to the time grid imposed by the raster plot 
measurement. In some sense this produces a stroboscopic dynamical system, with a frequency given by 
the time resolution 5a, producing a phenomenological representation of the underlying continuous time 
evolution. 

This has several advantages, (i) this simplifies the mathematical analysis of the dynamics avoiding 
the use of delta distributions, left-right limits, etc ... appearing in the continuous version; (ii) provided 
that mathematical results do not depend on the finite time discretization scale, one can take it arbitrary 
small; (iii) it enhances the role of symbolic coding and raster plots. 

Henceforth, from now on, we fix a positive time scale 5 > which can be mathematically arbitrary 
small, such that (i) a neuron can fire at most once between [t, t + 5[ (i.e. 5 « r, the refractory period); 
(ii) dt << 5, so that we can keep the continuous time evolution of membrane potentials (3), taking 
into account time scales smaller than 5, and integrating membrane potential dynamics on the intervals 
[t,t + 5[; (iii) the spike time is known within a precision 5. Therefore, the terminology, "neuron k fires 
at time t" has to be replaced by "neuron k fires between t and t + 5"; (iv) the update of conductances is 
made at times multiples 8 of 5. 

1.2.2 Raster plot. 

In this context, we introduce a notion of "raster plot" which is essentially the same as in biological 

measurements. A raster plot is a sequence u> d = {<*>(t)}J^, of vectors u>(t) = f [u)k{t)]^ =1 such that the 
entry Wfc(t) is 1 if neuron k fires between [t, t + S[ and is otherwise. Note however that for mathematical 

8 This could correspond to the following "experiment". Assume that we measure the spikes emitted by a set of in vitro 
neurons, and that we use this information to update the conductances of a model like (5), in order to see how this model 
"matches" the real neurons (see (Jolivet, Rauch, Lscher, & Gerstner, 2006) for a nice investigation in this spirit). Then, we 
would have to take into account that the information provided by the experimental raster plot is discrete, with a clock-based 
grid, even if the membrane potential evolves continuously. 



9 



reasons, explained later on, a raster plot corresponds to the list of firing states {u(t)}^ over an infinite 
time horizon, while on practical grounds one always considers bounded times. 

Now, for each k — 1...N, one can decompose the interval I = [V m i n , V max ] into 1 U Ti with 
1q = \V m in, = [9, V m ax\- If Vk £ %o neuron k is quiescent, otherwise it fires. This splitting induces 

a partition V of M, that we call the "natural partition". The elements of V have the following form. 
Call A = {0, 1}^. Let lj = [uk\k=i £ A. This is a TV dimensional vector with binary components 0, 1. 

We call such a vector a firing state. Then M = [J where: 

u>eA 

Mu = {VgM I V k e J W J. (9) 

Therefore the partition V corresponds to classifying the membrane potential vectors according to 
their firing state. Indeed, to each point V(i) of the trajectory V corresponds a firing state u>{t) whose 
components are given by: 

w fc (t) = Z[V k (t)], (10) 

where Z is defined by : 

Z{x) = X [x>B], (11) 

where \ is the indicator function that will later on allows us to include the firing condition in the evolution 
equation of the membrane potential (see (20)). On a more fundamental ground, the introduction of raster 
plots leads to a switch from the dynamical description of neurons, in terms of their membrane potential 
evolution, to a description in terms of spike trains where & provides a natural "neural code" . From the 
dynamical systems point of view, it introduces formally a symbolic coding and symbolic sequences are 
easier to handle than continuous variables, in many aspects such as the computation of topological or 
measure theoretic quantities like topological or Kolmogorov-Sinai entropy (Katok & Hasselblatt, 1998). A 
natural related question is whether there is a one-to-one correspondence between the membrane potential 
trajectory and the raster plot (see theorem 2). 

Note that in the deterministic models that we consider here, the evolution, including the firing times 
of the neurons and the raster plot, is entirely determined by the initial conditions. Therefore, there is no 
need to introduce an exogenous process (e.g. stochastic) for the generation of spikes (see (Destexhe & 
Contrcras, 2006) for a nice discussion on these aspects). 

Furthermore, this definition has a fundamental consequence In the present context, current and con- 
ductances at time t become functions of the raster plot up to time t. Indeed, we may write (7) in the 
form: 

Mj(t.w) 

gk (t,[Q] t )^Y.G% E ^{s-tf^ds (12) 

j n=l 

where [Q] t = {w(s)}* =0 is the raster plot up to time t and Mj(t,uo) is the number of spikes emitted by 
neuron j up to time t, in the raster plot u> (i.e. Mj{t,Cb) = Y^s^i^ji 8 ))- now ^ 1S a mmt ipl c °f 

Remark. In continuous time Integrate and Fire models u> can assume uncountably many values. 
This is because a neuron can fire at any time and because firing is instantaneous. Therefore, the same 
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property holds also if one considers sequences of firing states over a bounded time horizon. This is still 
the case even if one introduces a refractory period, because even if spikes produced by a given neuron are 
separated by a time slot larger or equal than the refractory period, the first spike can occur at any time 
(with an infinite precision). If, on the opposite, we discretize time with a time scale 5 small enough to 
ensure that each neuron can fire only once between t and t + S, ui, truncated at time T5 can take at most 
2 NT values. For these reasons, the "computational power" of Integrate and Fire models with continuous 
time is sometimes considered as infinitely larger than a system with clocked based discretization (Maass 
& Bishop, 2003). The question is however whether this computational power is something that real 
neurons have, or if we are dealing with a model-induced property. 

1.2.3 Integrate regime. 

For this regime, as we already mentioned, we keep the possibility to have a continuous time (dt « 5) 
evolution of membrane potential (3). This allows us to integrate V on time scales smaller than 5. But, 
since conductances and currents depends now on the raster plot cD, we may now write (3) in the form: 

dV k 

—— + g k (t, [w] t ) V k = ik{t, [to] t ); whenever V k < 8. (13) 

When neuron k does not fire between t,t + 5 one has, integrating the membrane potential on the 
interval t, t + 5 (see appendix) : 

V k (t + S) = lk (t, [u>) t ) V k (t) + J k (t, [u>] t ), (14) 

where: 

lk (t,[Q} t )^e-C 5 MQ] > )ds , (15) 

and where: 

rt+6 

J k (t,[u>] t ) = J ik(s,[Q] t )vk{s,t + S,[u)] t )ds, (16) 
is the integrated current with: 

v k {s,t + 5Mt) = e-^ +S 3k(s ' A ^ )ds '. (17) 

Remarks. 

1. In the sequel, we assume that the external current (see (8)) is time- constant. Further developments 
with a time dependent current, i.e. in the framework of an input-output computation (Bertschingcr 
& Natschlager, 2004), will be considered next. 

2. We note the following property, central in the subsequent developments. Since g k (t, [<D] t ) > 0, 

7fc(*,[w] t ) < l,Vt, Vw, Vfc. (18) 
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1.2.4 Firing regime. 



Let us now consider the case where neuron k fires between t and t + 8. In classical IF models this means 
that there is some ij^ S [t,£ + <J[ such that V k {t^) = 9. Then, one sets V k (t k n ^ + ) — V reset (instantaneous 
reset). This corresponds to Fig. lb. Doing this one makes some error compared to the real spike shape 
depicted in Fig. la. In our approach, one does not know exactly when firing occurs but we use the 
approximation that the spike is taken into account at time t + 8. This means that we integrate V k until 
t + 5 then reset it. In this scheme V k can be larger than 9 as well. This explains why Z(x) = x(x > 9). 
This procedures corresponds to Fig. lc (alternative I). One can also use a slightly different procedure. 
We reset the membrane potential at t + 8 but we add to its value the integrated current between [t, t + 8[. 
This corresponds to Fig. Id (alternative II). We have therefore three types of approximation for the real 
spike in Fig. la. Another one was proposed by (Hansel et al., 1998), using a linear interpolation scheme. 
Other schemes could be proposed as well. One expects them to be all equivalent when 8 — > 0. For finite 
8, the question whether the error induced by these approximations is crucial and discussed in section 2.6. 

In this paper we shall concentrate on alternative II though mathematical results can be extended to 
alternative I in a straightforward way. This corresponds to the initial choice of the Beslon-Mazet-Soula 
Model motivating the paper (Soula et al., 2006) and the present work. 

In this case, the reset corresponds to : 

V k (t)>0=>V k (t + S) = J k {t,[u] t ), (19) 

(recall that V rese t = 0). 

Integrate and Fire regime can now be included in a unique equation, using the function Z defined in 
(11): 

V k (t + 1) = 7k (t, [u>] t ) [1 - Z(V k (t))} V k (t) + J k (t, [Q] t ), (20) 
where we set 8=1 from now on. 

V(t) V(t) V(t) V(t) 




Figure 1: A. "Real spike" shape; the sampling window is represented at a scale corresponding to a 
"small" sampling rate to enhance the related bias. B. Spike shape for an integrate and fire model with 
instantaneous reset, the real shape is in blue. C. Spike shape when reset occurs at time t + 8 (alternative 

I) . D. Spike shape with reset at time t + 8 plus addition of the integrate current (green curve) (alternative 

II) . 
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1.3 Examples. 

1.3.1 The Beslon-Mazet-Soula Model (BMS). 

Consider the leaky integrate and fire model, where conductances are constant. Set W k j = G^E + (W k j — 
G^-E~) for excitatory (inhibitory) synapses. Then, replacing the a-profile by a Dirac distribution, (20) 
reduces to: 

N 

V k (t + l)= 7 V k (t) [l-Z(V k (t))]+Y / W kl Z(V J (t)) + l { ext) (21) 

This model has been proposed by G. Beslon, O. Mazet and H. Soula in (Soula et al., 2006). A 
mathematical analysis of its asymptotic dynamics has been done in (Cessac, 2008) and we extend these 
results to the more delicate case of conductance based IF models in the present paper. (Note that having 
constant conductances leads to a dynamics which is independent of the past firing times (raster plot). In 
fact, the dynamical system is essentially a cellular automaton but with a highly non trivial dynamics). 

1.3.2 Alpha-profile conductances. 

Consider now a conductance based model of form (3), leading to: 

7fe(t, [uj) t ) = Ke ^ *AZ^„ = i i t ^ o i ) } (22) 

where if is a constant: 

K = e~4; < 1. (23) 

while, using the form (6) for a gives: 

._,(») 



lk {t,[u] t ) = Ke l L J n (24) 

One has therefore to handle an exponential of an exponential. This example illustrates one of the main 
problem in IF models. IF models have been introduced to simplify neurons description and to simplify 
numerical calculations (compared e.g. with Hodgkin-Huxley's model (Hodgkin & Huxley, 1952)). Indeed, 
their structure allows one to write an explicit expression for the next firing times of each neurons, knowing 
the membrane potential value. However, in case of a exponential profile, there is no simple form for the 
integral and, even in the case of one neuron, one has to use approximations with T functions (Rudolph & 
Destexhe, 2006) which reduce consequently the interest of IF models and event based integration schemes. 

2 Theoretical analysis of the dynamics. 
2.1 The general picture. 

In this section we develop in words some important mathematical aspects of the dynamical system (20), 
mathematically proved in the sequel. 
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Singularity set. The first important property is that the dynamics (20) (and the dynamics of 
continuous time IF models as well) is not smooth, but has singularities, due to the sharp threshold 
definition in neurons firing. The singularity set is: 

S = {V e M\3i = 1 . . . N, such that V. t = 9} . 

This is the set of membrane potential vectors such that at least one of the neurons has a membrane 
potential exactly equal to the threshold 9 . This set has a simple structure: it is a finite union of TV — 1 
dimensional hypcrplanes. Although S is a "small" set both from the topological (non residual set) and 
metric (zero Lebesgue measure) point of view, it has an important effect on the dynamics. 

Local contraction. The other important aspect is that the dynamics is locally contracting, because 
"fk(t, p] t ) < 1 (see eq. (18)). This has the following consequence. Let us consider the trajectory of a 
point V 6 M and perturbations with an amplitude < e about V (this can be some fluctuation in the 
current, or some additional noise, but it can also be some error due to a numerical implementation). 
Equivalcntly, consider the evolution of the e-ball £>(V, e). If £>(V, e) n S = then we shall see in the 
next section that the image of B(V, e) is a ball with a smaller diameter. This means, that, under the 
condition S(V, e) PI S = 0, a perturbation is damped. Now, if the images of the ball under the dynamics 
never intersect S, any e-perturbation around V is exponentially damped and the perturbed trajectories 
about V become asymptotically indistinguishable from the trajectory of V. Actually, there is a more 
dramatic effect. If all neurons have fired after a finite time t then all perturbed trajectories collapse onto 
the trajectory of V after t + 1 iterations (see prop. 1 below). 

Initial conditions sensitivity. On the opposite, assume that there is a time, to, such that the image 
of the ball B(V, e) intersects S. By definition, this means that there exists a subset of neurons {i\, . . . , ik} 
and V 6 B(V,e), such that Z(Vi(to)) ^ Z(V[(to)), i <E {ii, ■ ■ ■ , ik}- For example, some neuron does 
not fire when not perturbed but the application of an e-perturbation induces it to fire (possibly with a 
membrane potential strictly above the threshold). This requires obviously this neuron to be close enough 
to the threshold. Clearly, the evolution of the unperturbed and perturbed trajectory may then become 
drastically different (see Fig. 2). Indeed, even if only one neuron is lead to fire when perturbed, it may 
induce other neurons to fire at the next time step, etc ... , inducing an avalanche phenomenon leading to 
unpredictability and initial condition sensitivity 10 . 

It is tempting to call this behavior "chaos" , but there is an important difference with the usual notion 
of chaos in differentiable systems. In the present case, due to the sharp condition defining the threshold, 
initial condition only occurs at sporadic instants, whenever some neuron is close enough to the threshold. 
Indeed, in certain periods of time the membrane potential typically is quite far below threshold, so that 
the neuron can fire only if it receives strong excitatory input over a short period of time. It shows then 
a behavior that is robust against fluctuations. On the other hand, when membrane potential is close to 
the threshold a small perturbation may induce drastic change in the evolution. 

Stability with respect to small perturbations. Therefore, depending on parameters such as 
the synaptic weights, the external current, it may happen that, in the stationary regime, the typical 
trajectories stay away from the singularity set, say within a distance larger than e > 0, which can be 

9 A sufficient condition for a neuron i to fire at time t is Vj(t) = 9 hence V(t) £ <S. But this is not a necessary condition. 
Indeed, as pointed in the footnote 5, there may exist discontinuous jumps in the dynamics, even if time is continuous, either 

i — — 

due to noise, or a profiles with jumps (e.g. a(t) = —e T , t > 0). Thus neuron i can fire with Vi(t) > 9 and V(t) f S. In 
the present case, this situation arises because time is discrete and one can have V(t — S) < 9 and V(t) > 9. This holds as 
well even if one uses numerical schemes using interpolations to locate more precisely the spike time (Hansel et al., 1998). 

10 This effect is well known in the context of synfire chains (Abeles, 1982, 1991; Abeles et al., 1993; Hertz, 1997) or 
self-organized criticality (Blanchard et al., 2000). 
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arbitrary small, (but positive). Thus, a small perturbation (smaller than e) does not produce any change 
in the evolution. At a computational level, this robustness leads to stable input-output transformations. 
In this case, as we shall see, the dynamics of (20) is asymptotically periodic (but there may exist a large 
number of possible orbits, with a large period). In this situation the system has a vanishing entropy 11 . 
This statement is made rigorous in theorem 1 below. 

On the other hand, if the distance between the set where the asymptotic dynamics lives 12 and the 
singularity set is arbitrary small then the dynamics exhibit initial conditions sensitivity, and chaos. Thus, 
a natural question is: is chaos a generic situation ? How does this depend on the parameters ? A related 
question is: how does the numerical errors induced by a time discretization scheme evolve under dynamics 
(Hansel et al., 1998) ? 

"Edge of chaos". It has been shown, in (Cessac, 2008) for the BMS model, that there is a sharp 
transition 13 from fixed point to complex dynamics, when crossing a critical manifold usually called the 
"edge of chaos" in the literature. While this notion is usually not sharply defined in the Neural Net- 
work literature, we shall give a mathematical definition which is moreover tractable numerically. Strictly 
speaking chaos only occurs on this manifold, but from a practical point of view, the dynamics is in- 
distinguishable 14 from chaos, close to this manifold. When the distance to the edge of chaos further 
increases the dynamics is periodic with typical periods compatible with simulation times. This manifold 
can be characterized in the case where the synaptic weights are independent, identically distributed with 

2 

a variance j^. 

In BMS model (e.g., time discretized gIF model with constant conductances) it can be proved that 
the chaotic situation is non generic (Cessac, 2008). We now develop the same lines of investigation 
and discuss how these result extend to the model (20). Especially, the "edge of chaos" is numerically 
computed and compared to the BMS situation. 

Let us now switch to the related mathematical results. Proofs are given in the appendix. 

2.2 Piecewise affine map. 

Let us first return to the notion of raster plot developed in section 1.2.2. At time t, the firing state 
u>(t) 6 A can take at most 2 N values. Thus, the list of firing states w(0), . . . ,u>(t) E A t+1 can take at 
most 2 w ( t+1 ) values. (In fact, as discussed below, only a subset of these possibilities is selected by the 
dynamics). This list is the raster plot up to time t and we have denoted by [<DL. Once the raster plot 
up to time t has been fixed the coefficients 7^ and the integrated currents Jk in (20) are determined. 
Fixing the raster plot up to time t amounts to construct branches for the discrete flow of the dynamics, 
corresponding to sub-domains of M constructed iteratively, via the natural partition (9), in the following 
way. 

Fix t > and [uj] t . Note: 

M [a]t = {V e M\V(s) e M uw> .=„...*} . 

11 More precisely the topological entropy (average bit rate production considered over an infinite time horizon) is zero 
but this implies that the Shannon entropy is also zero. 

12 Say the "attractor", though one must be cautious with this notion, as we shall see below. 

13 This transition is reminiscent of the one exhibited in (Keener et al., 1981) for an isolated neuron submitted to a periodic 
excitation, but the analysis in (Cessac, 2008) and the present analysis hold at the network level. 

14 Namely, though the dynamics is periodic, the periods are well beyond the times numerically accessible. 



15 



This is the set (possibly empty) of initial membrane potentials vectors V = V(0) whose firing pattern 
at time s is w(s), s = . . .t. Consequently, VV e , we have: 

t t t 

V k (t+l) = l[ lk ( S ,^} s )[l-u; k ( S )}V k (0) + Y / MnA^ n ) II 7*(«,N.)[l-WfcW] k = l...N, (25) 

s— n— s= n+1 

as easily found by recursion on (20). We used the convention n'it 7fc( s ; PL) [1 — w fc( s )] = 1- 
Then, define the map: 

F ^ "\ V - F^ +1 V = V(t + l) ' (26) 
with Vfe(t + 1) given by (25) and F% = Id. Note that F^ +1 is affine. Finally define: 

f F*+ 1 :A< - JW 

\ V6% - F^(V) ^ 

such that the restriction of F t+1 to the domain M-p] is precisely F^ +1 . This mapping is the flow of the 
model (20) where: 

V(t + 1) = F t+1 V, VeM. 

A central property of this map is that it is piecewise affine and it has at most 2 Nlyt+1 ^ branches 
F^ +1 parametrized by the legal sequences [Cj] t which parametrize the possible histories of the conduc- 
tance/currents up to time t. 

Let us give a bit of explanation of this construction. Take V = V(0) G Mu>(o)- This amounts 
to fixing the firing pattern at time with the relation ui k (0) = Z(V k (0)), k = 1...N. Therefore, 
Vfe(l) = 7fc(0,w(0)) [1 - Wfe(0)] Vfe(0) + J fe (0,u>(0)), where j k ,J k do not depend on V(0) but only on 
the spike state of neurons at time 0. Therefore, the mapping Fi : A^ W (o) — * -M such that F^.~V = 
7fc(0,u))[l — Wfc (0)] Vfc (0) + J k (0,u)),k = 1...N is affine (and continuous on the interior of Mu(o))- 
Since ^(0) is an hypercube, FiA^t^o) is a convex connected domain. This domain typically intersects 
several domains of the natural partition V . This corresponds to the following situation. Though the 
pattern of neuron firing at time is fixed as soon as V(0) £ M-u)(Q)i the list of neurons firing at the next 
time depends on the value of the membrane potentials V(0), and not only on the spiking pattern at time 
0. But, by definition, the domain: 

Mu>(i)u>(o) = Wp], = (Fi ) 1 M U (i) n Mu3(a) 

is such that VV(0) G ■M.u(i)u(o)> the spiking pattern at time is w(0) and it is w(l) at time 1. If 
the intersection is empty this simply means that one cannot find a membrane potential vector such that 
neurons fire according to the spiking pattern u>(0) at time t = then fire according to the spiking pattern 
w(l) at time t = 1. If the intersection is not empty we say that "the transition u>(0) — > w(l) is legal 15 " . 

Proceeding recursively as above one constructs a hierarchy of domains M-[ui\ t and maps F^ +1 . Inci- 
dentally, an equivalent definition of M[a>] t is: 

15 Conditions ensuring that a transition is legal depend on the parameters of the dynamical systems, such as the synaptic 
weights. 
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Figure 2: Schematic representation, for two neurons, of the natural partition V and the mapping discussed 

in the text. In this case, a firing state is a vector with components u) = (j*>%) labeling the partition 

elements. A set of initial conditions, say a small (L°°) ball in A4u>, is contracted by leak (neuron 1 in 
the example) and reset (neuron 2 in the example), but its image can intersect the singularity set. This 
generates several branches of trajectories. Note that we have given some width to the projection of the 
image of the ball on direction 2 in order to see it on the picture. But since neuron 2 fires the width is in 
fact 0. 



t 

M"] t = rK F *r lA w ( 28 ) 

s=0 

As stated before, is the set of membrane potential vectors V such that the firing patterns up 

to time t are u>(0), ■ ■ ■ , If this set is non empty we say that the sequence k>(0), ■ • ■ ,<*>(£) is legal. 

Though there are at most 2 jv ^ t+1 ' possible raster plots at time t the number of legal raster plots is 
typically smaller. This number can increase either exponentially with t or slower. We shall denote by 
£^ the set of all legal (infinite) raster plots (legal infinite sequences of firing states). Note that is a 
topological space for the product topology generated by cylinder sets (Katok & Hasselblatt, 1998). The 
set [uj] t of raster plots having the same first t + 1 firing patterns is a cylinder set. 

2.3 Phase space contraction. 

Now, we have the following: 

Proposition 1 For all Co £ S^,Vi > 0, the mapping V <E M.[q] — > F^ +1 (V) is affine, with a Jacobian 
matrix and an affine constant depending on t,[u)] t . Moreover, the Jacobian matrix is diagonal with 
eigenvalues 

t 

a k (t, u>) = JJ 7 fe(s, [£>U(1 " WkOO) < 1, k = 1 . . .N, 

s=0 
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Consequently, (V) is a contraction. 

Proof The proof results directly from the definition (26) and (25) with 7fc(s, [u>] s ) < l,Vs > (see (18) 
)• ■ 

Since the domains A4uj of the natural partition are convex and connected, and since F is affine on 
each domain (therefore continuous on its interior), there is a straightforward corollary: 

Corrollary 1 The domains M.[u\ are convex and connected. 

There is a more important, but still straightforward consequence: 

Corrollary 2 F* +1 is a non uniform contraction on M where the contraction rate in direction k is 

^El= logM^)],VVeA^ ]t . 

Then, we have the following : 
Proposition 2 Fix u> G E + . 

1. If 3t < oo, such that, Vfc = 1 . . . N, 3s = s(k) < t where u)k(s) = 1 then F~ + ] is a point. 
That is, all orbits born from the domain A4[q] converge to the same orbit in a finite time. 

2. If 3k 6 {1,...N} such that Vt > 0, Wfc(i) = then F^ +1 is contracting in direction k, with a 
Lyapunov exponent Afc(cD), such that: 

1 ' 1 * 

liminf — — - Vlog7 fc (s, < X k (w) < limsup — — - V log7 fc (s, [w] s ) < 

s— s— 

Proof Statement 1 holds since, under these hypothesis, all eigenvalues of F^ +1 arc zero. For 2, since 
DF^ +1 is diagonal, the Lyapunov exponent in direction k is defined by A& (w) = lim^oo S«=o logC' 7 *: (*> ^ 
whenever the limit exists (it exists almost surely for any F invariant measure from Birkoff theorem). ■ 

Remark. An alternative definition of Lyapunov exponent has been introduced by Coombes in 
(Coombes, 1999b) and (Coombes, 1999a), for Integrate and Fire neurons. His definition, based on 
ideas developed for impact oscillators (Muller, 1995), takes care of the discontinuity in the trajectories 
arising when crossing S. Unfortunately, his explicit computation at the network level (with continuous 
time dynamics), makes several implicit assumptions (see eq. (6) in (Coombes, 1999a)): (i) there is a finite 
number of spikes within a bounded time interval; (ii) the number of spikes that have been fired up to 
time t, Vt > 0, is the same for the mother trajectory and for a daughter trajectory, generated by a small 
perturbation of the mother trajectory at t — 0; (iii) call Xf, in Coombes' notations, the fc-th spike time 
for neuron i in the mother trajectory, and the fc-th spike time for neuron i in the daughter trajectory. 
Then T* = T^ + Sf, where 6"? is assumed to become arbitrary small, Vfc > 0, when the initial perturbation 
amplitude tends to zero. While assumption (i) can be easily fulfilled (e.g. by adding a refractory period) 
assumptions (ii) and (iii) are more delicate. 

Transposing this computation to the present analysis, this requires that both trajectories are never 
separated by the singularity set. A sufficient condition is that the mother trajectory stays sufficiently 
far from the singularity set. In this case the Lyapunov exponent defined by Coombes coincides with our 
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definition and it is negative. On the other hand, in the "chaotic" situation (see section 2.4.3), assumptions 
(ii) and (iii) can typically fail. For example, it is possible that neuron i stops firing after a certain time, 
in the daughter trajectory, while it was firing in the mother trajectory, and this can happen even if 
the perturbation is arbitrary small. This essentially means that the explicit formula for the Lyapunov 
exponent proposed in (Coombes, 1999a) cannot be applied as well in the "chaotic" regime. 

2.4 Asymptotic dynamics. 

2.4.1 Attracting set A and w-limit set. 

The main notion that we shall be interested in from now on concerns the invariant set where the asymp- 
totic dynamics lives. 

Definition 1 (From (Guckenheimer & Holmes, 1983; Katok & Hasselblatt, 1998)) A point y G Ai is 
called an u-limit point for a point x G M. if there exists a sequence of times {tk\X=o> suca ^at x(tk) — > y 
as tk — > +oo. The uj-limit set of x, u)(x), is the set of all uo-limit points of x. The iv-limit set of M, 
denoted by O, is the set Q = \J xeM uj(x). 

Equivalently, since M is closed and invariant, we have 12 = Ht^o F'(A4). 

Basically, ft is the union of attractors. But for technical reasons, related to the case considered in 
section 2.4.3, it is more convenient to use the notion of w-limit set. 

2.4.2 A theorem about the structure of 0. 

Theorem 1 Assume that Be > and 3T < oo such that, VV G M, Vi G {1 

1. Either 3t < T such that Vi{t) > 9; 

2. Or 3t Q = t (V, e) such that Vi > t , Vi(t) < 9 - e 
Then, ft is composed by finitely many periodic orbits with a period < T. 

The proof is given in the appendix A. 2. 

Note that conditions (1) and (2) are not disjoint. The meaning of these conditions is the following. (1) 
corresponds to imposing that a neuron has fired after a finite time T (uniformly bounded, i.e. independent 
of V and i). (2) amounts to requiring that after a certain time to, the membrane potential stays below the 
threshold value and it cannot accumulate on 9. We essentially want to avoid a situation when a neuron 
can fire for the first time after an unbounded time (see section 2.4.3 for a discussion of this case). Thus 
assumptions (1),(2) look quite reasonable. Under these assumptions the asymptotic dynamics is periodic 
and one can predict the evolution after observing the system on a finite time horizon T, whatever the 
initial condition. Note however that T can be quite a bit large. 

There is a remarkable corollary result, somehow hidden in the proof given in the appendix. The 
neurons that do not fire after a finite time are still driven by the dynamics of firing neurons. It results 
that, in the asymptotic regime, non firing neurons have a membrane potential which oscillates below the 
threshold. This exactly corresponds to what people call "sub-threshold oscillations" . In particular, there 
are times where those membrane potentials are very close to the threshold, and a small perturbation 
can completely changes further evolution. This issue is developed in the next section. Possible biological 
interpretations are presented in the discussion section. 



...N}, 
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2.4.3 Ghost orbits. 



The advantage of the previous theorem is that we define conditions where one can relatively easily controls 
dynamics. However, what happens if if we consider the complementary situation corresponding to the 
following definition ? 

Definition 2 An orbit V is a ghost orbit if 3i such that: 

(i) Vi > 0,Vi(t) < 9 

and : 

(ii) limsupl^(i) = 9 

t — >+oc 

Namely there exists at least one initial condition V and one neuron i such that one cannot uniformly 
bound the first time of firing, and Vi{t) approaches arbitrary close the threshold. In other words sub- 
threshold oscillations drive the neuron "dangerously close" to the threshold, though we are not able to 
predict when the neuron will fire. This may look a "strange" case from a practical point of view, but 
it has deep implications. This indeed means that we can observe the dynamics on arbitrary long times 
without being able to predict what will happen later on, because when this neuron eventually fire, it may 
drastically change the evolution. This case is exactly related to the chaotic or unpredictable regime of 
IF models. 

One may wonder whether the existence of ghost orbits is "generic" . To reply to this question one 
has first to give a definition of genericity. In the present context, it is natural to consider the dynamical 
system describing the time evolution of our neural network as a point in a space H of parameters. These 
parameters can be e.g., the synaptic weights, or parameters fixing the time scales, the reversal potentials, 
the external currents, etc... Varying these parameters (i.e., moving the point representing our dynamical 
system in H) can have two possible effects. Either there is no qualitative change in the dynamics and 
observable quantities such as e.g., firing rates, average inter-spikes interval, etc, are varying continuously. 
Or, a sharp change (bifurcation) occurs. This corresponds to the crossing of a critical or bifurcation 
manifold in Tt. Now, a behavior is generic if it is "typical" . On mathematical grounds this can have two 
meanings. Either this behavior is obtained, with a positive probability, when drawing the parameters 
(the corresponding point in H) at random with some natural probability (e.g., Gaussian). In this case 
one speaks of "metric genericity" . Or, this behavior holds in a dense subset of TL. One speaks then of 
"topological genericity" . The two notions usually do not coincide. 

In the BMS model, ghost orbits arc non generic in both senses (Cessac, 2008). The proof does not 
extend to more general models such as (20) because it heavily uses the fact that the synaptic current 
takes only finitely many values in the BMS model. As soon as we introduce a dependence in Cj this is 
not the case anymore. We don't know yet how to extend this proof. 

2.5 Edge of chaos. 

On practical grounds ghost orbits involve a notion of limit t — > +oo which has no empirical meaning. 
Therefore the right question is: are there situations where a neuron can fire for the first time after a time 
which is well beyond the observation time ? One way to analyze this effect is to consider how close the 
neurons are to the threshold in their evolution. On mathematical grounds this is given by the distance 
from the singularity set to the w-limit set: 
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d(Q,S)= inf inf min IVM - 0\. (29) 
vent>oi=i...Af 

The advantage of this definition, is that it can easily be adapted to the plausible case where observa- 
tion time is bounded (see section 3). 

Now, the following theorem holds. 
Theorem 2 . 

1. If d(Cl,S) > then il is composed by finitely many periodic orbits with a finite period. 

2. There is a one-to-one correspondence between a trajectory on Q and its raster plot. 

The proof is exactly the same as in (Cessac, 2008) so we don't reproduce it here. It uses the fact that 
if d(Cl, S) > then there is a finite time T, depending on d(Cl, S) and diverging as d(Cl, S) —> 0, such that 
F T has a Markov partition (constituted by local stable manifolds since dynamics is contracting) where 
the elements of the partition are the domains M[&] . Note however that d(Q,S) > is a sufficient but 
not a necessary condition to have a periodic dynamics. In particular, according to theorem 1 one can 
have d(Q, S) = and still have a periodic dynamics, if at some finite time t, for some neuron i, Vi(t) = 9. 
This strict equality is however not structurally stable, since a slight change e.g. in the parameters will 
remove it. The main role of the condition d(Cl, S) > is therefore to avoid situations where the membrane 
potential of some neuron accumulates on 9 from below (ghost orbits). See the discussion section for a 
possible biological interpretation on this. 

But d(Cl,S) plays another important role concerning the effect of perturbations on the dynamics. 
Indeed, if d(fl,S) > then the contraction property (corollary 2) implies that any perturbation smaller 
than G?(f2,5) will be damped by dynamics. In particular, making a small error on the spike time, or any 
other type of error leading to an indeterminacy of V smaller than d(Cl, S) will be harmless. 

Let us now discuss the second item of theorem 2. It expresses that the raster plot is a symbolic coding 
for the membrane potential trajectory. In other words there is no loss of information on the dynamics 
when switching from the membrane potential description to the raster plot description. This is not true 
anymore if d(Cl,S) = 0. 

The first item tells us that dynamics is periodic, but period can be arbitrary long. Indeed, following 
(Cessac, 2008) an estimate for an upper bound on the orbits period is given by: 

n M ^ 2 N i° g «-y» (30) 

where < 7 > denotes the value of 7 averaged over time and initial conditions 16 (see appendix for details). 
Though this is only an upper bound this suggests that periods diverge when d(Q,S) — > 0. In BMS 
model, this is consistent with the fact that when d(Q, S) is close to dynamics "looks chaotic" . There- 
fore, d(Cl,S) is what a physicist could call an "order parameter", quantifying somehow the dynamics 
complexity. The distance d(fl,S) can be numerically estimated as done in (33) and (34), section 3. 

Before, we need the following list of (operational) definitions. 

Definition3 (Edge of chaos) 

The edge of chaos is the set of points £ in the parameter space TL where d(Q,S) = 0. 



3 Note that the system is not uniquely ergodic (see (Katok & Hasselblatt, 1998) for a definition of unique ergodicity). 



21 



The topological structure of £ can be quite complicated as we checked in the simplest examples (e.g., 
the BMS model with Laplacian coupling) and suggested by the papers (Bressloff & Coombes, 2000b, 
2000a) (see discussion). There are good reasons to believe that this set coincides with the set of points 
where the entropy is positive (see (Kruglikov & Rypdal, 2006a, 2006b) and discussion below). The set 
of points where the entropy is positive can have a fractal structure even in the simplest examples of one 
dimensional maps (Mackay & Tresser, 1986; Gambaudo & Tresser, 1988). Therefore, there is no hope to 
characterize £ rigorously in a next future. Instead, we use below a numerical characterization. 

The edge of chaos is a non generic set in the BMS model, and the same could hold as well in model 
(20). Nevertheless, it has a strong influence on the dynamics, since crossing it leads to drastic dynamical 
changes. Moreover, close to £ dynamics can be operationally indistinguishable from chaos. More precisely, 
let us now propose another definition. 

Definition 4 (Effective entropy) 

Fix T a called "the observational time". This is the largest accessible duration of an experiment. Call 
n(t) the number of (legal) truncated raster plots up to time t. Then, the effective entropy is; 

= i_logn(T ) (31) 

J- o 

Note that in the cases where raster plots provide a symbolic coding for the dynamics then liniT ^oo h^ e ^^ = 
n (top) ^ ^ ne topological entropy. 

On practical grounds, this definition corresponds to the following notion. The larger the effective 
entropy, the more the system is able to produce distinct neural codes. This provides one way to measure 
the "complexity" of the dynamics. On more "neuronal" grounds this quantity measures the variability 
in the dynamical response of the neuronal network to a given stimulus (external current) or its ability 
to produce distinct "functions" (a function being the response to a stimulus in terms of a spikes train) . 
Thinking of learning mechanisms (e.g., Hcbbian) and synaptic plasticity (LTD,LTP,STDP) one may 
expect to having the largest learning capacities when this entropy is large. This aspect will be developed 
in a separated paper. (For the effect of Hebbian learning and entropy reduction in firing rate neural 
networks see (Siri, Berry, Cessac, Delord, & Quoy, 2008)). 

Finally, a positive effective entropy means that the system essentially behaves like a chaotic system 
during the time of the experiment. Indeed, the entropy is closely related to the distance d(Q,S), since, 
from (30), a rough estimate/bound of h^ e ^^ is easily obtained from (30), (31): 

The notion of effective entropy provides some notion of "width" to the edge of chaos £ . For a fixed 
T Q the system behaves chaotically in a thick region £ To D £ in H such that h<- ef ^ > 0. And from (32) 
one expects that this entropy gets larger when d(Sl, S) gets smaller. 

2.6 Effects of time discretization. 

Under the light of the previous results, let us reconsider the approximation where spikes are taken into 
account at multiple of the characteristic time scale S, for the conductances update. Doing this, we 
make some error in the computation of the membrane potential, compared to the value obtained when 
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using the "exact" spike time value. Now, the question is whether this error will be amplified by the 
dynamics, or damped. As we saw, dynamics (20) is contracting but the effect of a small error can have 
dramatic consequences when approaching the singularity set. The distance d(fl,S) provides a criterion 
to define a "safe" region where a small error of amplitude e > in the membrane potential value is 
harmless, basically, if e < d(Q, S). On the other hand, if we are in a region of the parameters space where 
d(fl,S) — then a slight perturbation has an incidence on the further evolution. Since 5 can be arbitrary 
small in our theorems we have a good control on the dynamics of the continuous time IF models except 
at the edge of chaos where d(Q, S) = 0. This enhances the question of mathematically characterizing this 
region in the parameter space H. Note indeed that numerical investigations are of little help here since 
we are looking for a parameter region where the distance d(fl,S) defined as an asymptotic limit, has to 
be exactly zero. The problem is that even sophisticated schemes (e.g., event based) are also submitted 
to round off errors. Therefore, as a conclusion, it might well be that all numerical schemes designed to 
approximate continuous time Integrate and Fire models display trajectory sensitivity to spike time errors, 
when approaching d(Q,S) = 0. 

3 A numerical characterization of the "edge of chaos". 
3.1 A "coarse-grained" characterization. 

As mentioned in the previous section an exact analytic computation of the edge of chaos in the general 
case seems to be out of reach. However, a "coarse grained" characterization can be performed at the 
numerical level and possibly some analytic approximation could be obtained. For this, we choose the 
synaptic weights (resp. the synaptic conductances) (and/or the external currents) randomly, with some 
probability Vw (Pnext)), where W is the matrix of synaptic weights (Wij — E ± Gfj) and \( ext ) the vector of 
external currents (recall that external currents are time constant in this paper). A natural starting point 
is the use of Gaussian independent, identically distributed variables, where one varies the parameters, 
mean and variance, defining the probability definition (we call them statistical parameters, see (Cessac & 
Samuelides, 2007; Samuelides & Cessac, 2007) for further developments on this approach). Doing these, 
one performs a kind of fuzzy sampling of the parameters space, and one somehow expects the behavior 
observed for a given value of the statistical parameters to be characteristic of the region of W, i [ext) that 
the probabilities Vw, V-^t) weight (more precisely, one expects to observe a "prevalent" behavior in the 
sense of (Hunt, Sauer, & Yorke, 1992)). 

The idea is then to estimate numerically d(Cl, S) in order to characterize how it varies when changing 
the statistical parameters. As an example, in the present paper, we select conductances (resp. synaptic 
weights) randomly with a Gaussian probability with a fixed mean and a variance jj-, and we study the 
behavior of d(Q, S) when a is varying. Note that the neural network is almost surely fully connected. We 
compute numerically an approximation of the distance d(fl,S), where we fix a transient time T r and an 
observation time T and average over several initial conditions V^™) , n = 1 . . . Nqi for a given sample of 
synaptic weights. Then we perform an average over several synaptic weights samples W m , m = 1 . . . Nyy. 
In a more compact form, we compute: 




(33) 



where: 
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dt xp) = min min min \v} n) (T r + t) - 9\ . (34) 

In this way, we obtain a rough and coarse grained location of the edge of chaos where the distance 
d(n, S) vanishes. 

We have performed the following experiments with two control parameters. 

• Variance of random synaptic weights. We randomly select the synaptic strength which mod- 
ulates the synaptic conductance using a Gaussian distribution so that 80% of the synapses are 
excitatory and 20% inhibitory. The average standard-variation a is varied. The value a = 0.5 
corresponds, in our units, to what is chosen in the literature when considering the biological dis- 
persion in the cortex (e.g., (Rudolph & Destexhe, 2006)). Note however that, at the present stage 
of investigation, Dale's principle is not taken into account. 

• Membrane leak time-constant. As an additional control parameter we vary the membrane leak 
around the usual tl = 1 ■ • ■ 20ms values. This choice is two-fold. The value tl — 20ms corresponds 
to in-vitro measurement, while tl — > 1ms allows to represent in-vivo conditions in the cortex. On 
the other hand, it acts directly on the average contraction rate < 7 > which is a natural control 
parameter. 

Each simulation randomly selects the initial potential values in a — 70...30mT^ range. For each 
condition the simulation is run for Nqi = 100 initial conditions and Ny\> = 10 random selection of the 
synaptic weights. With a sampling period of 0.1ms, the network is run during T r — 1000 steps in order 
to "skip" the transients 17 and then observed during T Q = 1000 steps. In order to explore the role of 
history dependent conductances on the dynamics we considered different models from the biologically 
plausible IF model to BMS model. More precisely we explore four modeling variants: 

1. Model I, defined in (20). 

2. Model II (20) with a fixed 7. The evolution equation of membrane potentials is thus given by: 



V k (t + 1) = (7) V k (t) [1 - Z(V k {t))] + J k (t, [£j] t ) 

where the average (7) is the value observed during the numerical simulation of model I. Note that (7) 
depends on the parameters a, tl • The goal of this simulation is to check the role of the fluctuations 
of j(t, [u)] t ), controlling the instantaneous contraction rate, compared to a mean-field model where 
7(t, [u)] t ) is replaced by its average. This corresponds to what is called "current based" synapses 
instead of "conductance based" synapses in the literature, (see e.g. (Brette et al., 2007)). 

3. Model III (20) approximation with a fixed 7 and simplified synapses. The evolution equation of 
membrane potentials is given by: 



V k (t + 1) = (7) V k (t) [1 - Z(V k (t))} Z{V 3 {t 6+)) + 25" £ Gr. Z(V 3 (t <T)) 

3 3 

17 Note the transients depend on the parameters and on the distance to the singularity set too. In particular, one can 
have transients that are well beyond the current capacity of existing computers. Therefore, our procedure gives a rough 
localization of the edge of chaos. Analytic computation would give a more precise localization. 
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In addition to the previous simplification, we consider so-called "current jump" synapses where 
the synaptic input simply corresponds to an impulse, added to the membrane potential equation. 
Here the magnitude of the impulse and its delay 5~ — 2ms and S + — 10ms in order to keep both 
characteristics as closed as possible to the previous condition. 

4. Model IV (20) with a fixed 7 and instantaneous simplified synapses. The evolution equation of 
membrane potentials is given by: 

V k (t+l) =< 7 > V k (t) [1 - Z{V k (t))] +Y J W l] Z{V ] (t))) 

3 

where in addition to the previous simplification, the delay has been suppressed and where Wij = 
E^Gfy This last condition strictly correspond to the original BMS model (Soula et al., 2006). 

The results are given below. We have first represented the average value (7) for the model I in 
the range a G [0.01,1], T£ 6 [10,40]ms (see Fig. 3). The quantity related to the contraction rate, is 
remarkably constant (with small variations within the range [0.965,0.995]). 



99 

T *I_ 0.9B5 - 




Figure 3: Average value of 7 for model I a G [0.01, 1], T£ G [10,40]ms. The profile is very similar for 
other models. 

Then, we have considered the average value of the current J k {t, [u)] t ), averaged over time, initial 
conditions and neurons and denoted by / to alleviate the notations (Fig. 4), the logarithm of the 
distance d(£l,S) (Fig. 5), and the average Inter Spike Interval (ISI- Fig. 6), for the four models. The 
main observations are the following. Average current and Inter Spike Intervals have essentially the same 
form for all models. This means that these quantities are not really relevant if one wants to discriminate 
the various models in their dynamical complexity. 
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Figure 4: Average current I for the models I (top left)-II (top right)- III (bottom left)- IV (bottom right) 
with a € [0.01, 1], t l E [10, 40]ms. 



The observation of the distance d(£l,S) is quite more interesting. First, in the four models, the 
distance becomes very small when crossing some "critical region" in the plane tx,7. This region has 
a regular structure for the BMS model, but its structure seems more complex for (20). Note however 
that the numerical investigations used here do not allow us to really conclude on this point. The most 
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log[d(A,S)] i „ 5 log[d(A.S)] j . 5 




Figure 5: Average value of log(<i(fi, 5)) for the models I (top left)-II (top right)- III (bottom left)- IV 
(bottom right) with a 6 [0.01, 1], t l £ [10,40]ras. 



remarkable fact is that, in models III and IV, the distance increases when a increases beyond this region, 
while it does not in models I and II. This corresponds to the following observation. When the d(il,S) is 
small, one observes a complex dynamics with no apparent period. One naturally concludes to a chaotic 
regime. As we saw, strictly speaking it is in fact periodic but since periods are well beyond observable 
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Figure 6: Average value of the Inter Spike Interval for the models I (top left)-II (top right)- III (bottom 
left)- IV (bottom right), with a e [0.01, 1], t l e [10,40]ms. 



times, the situation is virtually chaotic 18 . When the distance increases, the orbits period decreases. 

18 Moreover, it is likely that the phase space structure has some analogies with spin-glasses. For example, if 7 = the 
dynamics is essentially equivalent to the Kauffman's cellular automaton (Kauffman, 1969). It has been shown by Derrida 
and coworkers (Derrida & Flyvbjerg, 1986; Derrida & Pomeau, 1986) that the Kauffman's model has a structure similar 
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Figure 7: Average value of log(d(f2, <S)) for the model I. (left) a E [1,10], r e [10,40] ms. We present 
here a projection in the plane cr, log d(£l, S), and the vertical bars correspond to the variations with tl- 
It allows us to verify the stability of the previous result for higher variability of the synaptic weights, 
(middle) t£ G [1, 1 . . . 2]ms below the usual 20ms value, a S [1, 10]. Such range corresponds to cortical 
neurons in high-conductance state. It allows to check the behavior of d(Q, S) in this case, (right) Sampling 
period of lms, in order to verify the robustness of the numerical results with respect to the sampling 
rate. 



to the Sherrington-Kirckpatrick spin-glass model(Mezard, Parisi, & Virasoro, 1987). The situation is even more complex 
when 7^0. It is likely that we have in fact a situation very similar to discrete time neural networks with firing rates where 
a similar analogy has been exhibited (Cessac, 1994, 1995). 29 



Therefore, there is a range of cr values where period become smaller than observational time and one 
concludes that dynamics is periodic. 

The situation is different for model I, II since the distance does not apparently increases with a. This 
suggests that introducing conductance based synapses and currents enhances considerably the width of 
the edge of chaos. On practical grounds, this means that model I, II have the capacity to display a 
very large number of distinct codes for wide choices of parameters. This is somewhat expected since 
the opposite conclusion would mean that introducing spike dependent conductances and current does 
not increases the complexity and information capacity of the system. But it is one thing to guess some 
behavior and another thing to measure it. Our investigations on the distance d(tt,S), a concept based 
on the previous mathematical analysis, makes a step forward in this direction. 

One step further, we have represented examples of raster plots in Fig. 8 and 9 for models I and IV. 
The figure 8 essentially illustrates the discussion above on the relation between the distance d(Q,S) and 
the dynamics; for a — 0.05, where d(fl,S) is "large", and dynamics is periodic; and for a = 0.4, where 
d(fl,S) is small, and dynamics looks more chaotic, for the two models. The difference between the two 
models becomes more accurate as a increases. Fig. 9 represents raster plots for models I, IV, with a = 10, 
where we study the effect of a small amount of noise, of amplitude 10~ 4 x 6 in the external current. This 
has no effect on model IV while it changes slightly the raster plot for model I, as expected. There is 
another remarkable difference. The code is sparser for model I than for model IV. This suggests that 
model I is in some sense optimal with respect to coding since it is able to detect very small changes in 
an input but the changes is not drastic and the neural code remains very sparse. 

4 Discussion 

We have thus an operational definition for the "edge of chaos" where an "order parameter" , the distance 
of orbits to the singularity has been defined. This parameter has a deep meaning. It controls how much 
the system is sensitive to perturbations. Such perturbations can be noise, but they can also be a small 
variation in the external current, corresponding e.g. to an input. If the amplitude of this perturbation is 
smaller than d(il, S) then it has no effect on the long term dynamics, and the neural code (raster plot) 
is unchanged. On the other hand, when the distance is small, even a tiny perturbation has a dramatic 
effect on the raster plot: the system produces a different code. As a corollary, the effective entropy is 
maximal when the distance is minimal. On practical ground, having a positive distance with a large 
effective entropy corresponds to situations where the system is able to produce a large number of distinct 
codes within the observational time, while this code is nevertheless robust to small perturbations of the 
input. Thus, we have a good compromise between the variability of the responses to distinct inputs and 
robustness of the code when an input is subject to small variations. 

Several questions are now open. A first one concerns the way how we measured this distance. We 
used a random sampling with independent synaptic weights. But these weights arc, in reality, highly 
correlated, via synaptic plasticity mechanism. What is the effect of e.g. STPD or Hebbian learning 
on the effective entropy is a perspective for a future work. Recent results in (Soula, 2005) and (Siri 
et al., 2008; Siri, Berry, Cessac, Dclord, & Quoy, 2007) suggest that synaptic plasticity reduces the 
entropy by diminishing the variability of raster plots and increasing the robustness of the response to 
an input. Some general (variational) mechanism could be at work here. This aspect is under investigation. 

Another important issue is the effect of noise. It is usual in neural network modeling to add Brownian 
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a = 0.05 



a = 0.40 



Figure 8: Examples of raster plots for the conductance based model (Model I, top row) and the leaky 
integrate and fire model (Model IV, bottom row). A time window of 100 samples is shown in each case. 
The control parameter is tl — 20ms. As visible in Fig. 5, a — 0.05 corresponds to a small order dynamics 
where the periodic behavior is clearly visible, and a = 0.40 to the "edge of chaos" . One blob width is 
lmsec 



noise to the deterministic dynamics. This noise accounts for different effects such as the diffusion of 
neurotransmitters involved in the synaptic transmission, the degrees of freedom neglected by the model, 
external perturbations, etc ... Though it is not evident that the "real noise" is Brownian, using this kind 
of perturbations has the advantage of providing a tractable model where standard theorems in the theory 
of stochastic processes (Touboul & Faugeras, 2007) or methods in non equilibrium statistical physics (e.g. 
Fokkcr-Planck equations (Brunei & Hakim, 1999)) can be applied. 

Though we do not treat explicitly this case in the present work, the formalism has been designed to 
handle noise effects as well. As a matter of fact, the effect of Brownian noise on the dynamics of our 
model can be analyzed with standard techniques in probability theory and stochastic perturbations of dy- 
namical systems (Freidlin & Wentzell, 1998). In particular, the probability distribution of the membrane 
potential trajectory can be obtained by using a discrete time version of Girsanov theorem (Samuelides & 
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without noise 



with noise 



Figure 9: Raster plots for models I (upper row) and IV (lower row), with a = 10.00 and the same 
condition as in Fig. 8. First column: model I and model without noise. Second column: same realization 
of synaptic weights and same initial conditions but with a small amount of noise in the external current. 
The noise is added to the membrane potential and its magnitude is very small (10~ 4 x 9). One blob 
width is 1msec. 



Cessac, 2007). Noise have several effects. Firstly, the stochastic trajectories stay around the unperturbed 
orbits until they jump to another attraction basin, the characteristic time depending on the noise inten- 
sity ("Arrhcnius law"). This has the effect of rendering the dynamics uniquely ergodic, which somehow 
simplifies the statistical analysis. The effect of noise will be essentially prominent in the region where 
d(Q,S) is small, leading to an effective initial condition sensitivity and an effective positive Lyapunov 
exponent, that could be computed using mean-field approaches (Cessac, 1995). It is possible to estimate 
the probability that a trajectory approaches the singularity set S within a finite time T and a distance 
d by using Freidlin-Wentzell estimates (Freidlin & Wentzell, 1998). One can also construct a Markov 
chain for the transition between the attraction basin of the periodic orbits of the unperturbed dynamics. 
The overall picture could be very similar (at least for BMS model) to what happens when stochastically 
perturbing Kauffman's model (Golinelli & Derrida, 1989), with possibly a phase space structure rcm- 
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iniscent of spin-glasses (where noise plays the role of the temperature). This study is under investigations. 

Yet another important issue relates to the fact that spikes can also be lost. This aspect is not yet 
taken into account in the present formalism, but annihilation of spikes is a future issue to address. 

A final issue is the relation of this work with possible biological observations. We would like in 
particular to come back to the abstract notion of ghost orbit. As said in the text, this notion corresponds 
to situation where the membrane potential of some "vicious" neuron fluctuates below the threshold, and 
approaches it arbitrary close, with no possible anticipation of its first firing time. This leads to an effective 
unpredictability in the network evolution, since when this neuron eventually fire, it may drastically change 
the dynamics of the other neurons, and therefore the observation of the past evolution does not allow 
one to anticipate what will be the future. In some sense, the system is in sort of a metastable state but 
it is not in a stationary state. 

Now, the biological intuition tends to consider that a neuron cannot suddenly fire after a very long 
time, unless its input changes. This suggests therefore that "vicious" neurons are biologically implausible. 
However, this argument, to be correct, must precisely define what is a "very long time". In fact, one has 
to compare the time scale of the experiment to the characteristic time where the vicious neurons will 
eventually fire. Note also that since only a very small portion of neurons can be observed e.g. in a given 
cortex area, some "vicious" neurons could be present (without being observed since not firing) , with the 
important consequence discussed in this paper. The observation of "temporarily silent" neurons which 
firing induces a large dynamic change would be an interesting issue in this context. 

As a final remark we would like to point out the remarkable work of Latham and collaborators 
discussing the effects induced by the addition or removal of a single spike in a raster plot. A central 
question is whether this "perturbation" (which is not necessarily "weak") will have a dramatic effect 
on the further evolution (see (Latham, Roth, Hausser, & London, 2006) and the talk of P. Latham 
available on line at http : / / www.archive.org/details/ Redwood JO enter _2006_09_25_Laift,am) . Especially 
the questions and discussions formulated during the talk of P. Latham are particularly salient in view of 
the present work. As an additional remark note that a perturbation may have an effect on trajectories 
but not on the statistics build on these trajectories (e.g. frequency rates) (Cessac & Sepulchre, 2006). 

A Appendix. 

A.l Computation of V k (t + 5) 
Fix t u t 2 e [t,t + 6[. Set: 

,, , r--, A - f' 2 g k (s,[^} 3 )ds - J* 2 g k (s,[u] t )ds 

v k {t 1 ,t 2 , H*) = e 1 = e 1 i ( 35 ) 

where the last equality holds from our assumption that spikes are taken into account at times multiples 
of S; therefore [u)] s — [u] t ,s£ [t,t + S[. 
We have: 

v k {t\,ti, = 1, 

Vk(ti,t 2 , [u>] t ) = ffe(*i, *i, [Lj] t )v k (t' 1 ,t 2 , [u)] t ), (36) 

for t[ e [t,t + 8[. Moreover: 
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dv k (t 1 ,t 2l [Q] t ) 



: g k (t lt [Q] t )v k (h,t2, [w] t ), 



This leads to: 



-7- {Vk(tl,t 2 , [w] t )Vfe(*l)) = v k (ti,t 2 , [u>] t ) 
at i 



^+g k (t ll [u] t )V k (t 1 ) 



v>k{ti,t 2 , [uj} t )ik(ti, [Q] t ). 



If neuron k does not fire between t and t+S we have, integrating the previous equation for t\ G [t, t + <5[ 
and setting i 2 = t + S : 



rt+6 

V k {t + S) =v k (t,t + 5,[u)} t )V k {t) + J i k {s,[u>] t )v k (s,t + 5,[u>] t )ds. (37) 



f-t+S 

It 

A. 2 Proof of theorem 1. 

The proof uses the following lemma. 

Lemma 1 Fix T a subset of {1 . . . TV} and let T be the complementary set of T ' . Call 



>,T,e = VeM 



(i) Vi G Tj 3t< T, such that V l {t) > 9 

(ii) Vj G 3t = t (V,j) < oo, such that Vi > t , Vj(t) <6-t 

then Q(rV,T,e)> the oj -limit set o/rV,T,e; is composed by finitely many periodic orbits with a period < T. 
Proof of theorem 1 

Note that there are finitely many subsets T of {1 . . . TV}. Note also that IY^e C IV ,t+i,e and that 
IV, T,e C Tjr T f , whenever e' < e. We have therefore: 

M c U U U T w = U rj7 .+°°.°- 

F T>0e>0 F 

But, under hypothesis (1) and (2) of theorem 1, there exists e > 0,T < oo such that M. — {JjrFf.T,e 
where the union on T is finite. Since ¥{M) G {Jjr F(IV,T,e), ft c (J^ n(Tjr T e y Under lemma 1 ft is 
therefore a subset of a finite union of sets containing finitely many periodic orbits with a period < T. ■ 

Proof of lemma 1 Call Hjr (resp. ILjr) the projection onto the subspace generated by the basis vectors 
ej, i G T (resp. e,, j G T) and set V> = n>V (V> = ILpV), = n>F (F^- = II^F). Since each 
neuron j G T is such that (25): 

t-i t-i 

= II 7fc(*,[w] a )(l- £ ,; fc (*))<(9-e ) (38) 

n— 1 s=n+l 

for i sufficiently large, (larger than the last (finite) firing time tj), these neurons do not act on the 
other neurons and their membrane potential is only a function of the synaptic current generated by the 
neurons G T . Thus, the asymptotic dynamics is generated by the neurons i G T. Then, VV G Cl(T^- t T t e), 
Vjr(t + 1) = F^[V^(t)] and V^(t + 1) = F^[V>(i)]- O ne can therefore focus the analysis of the u> limit 
set to its projection Q;F(IV,T,e) = ^-^{^^,T,e) (and infer the dynamics of the neurons j G T via (38)). 
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Construct now the partition P^ T \ with convex elements given by , where T is the same as in the 

definition of IV^e- By construction, F T+1 is continuous on each element P^ and fixing M[cj\ t amounts 

to fix the affinity constant of F T+1 . By definition of T, Z)F^- +1 | V , the derivative of F^- +1 at V, has all 

its eigenvalues equal to whenever V G f^IV^e) (prop. 1). Therefore F^- +1 [A^j n £V(IV,7\e)] i s a 
point. Since 

f? +1 (m n Mi>,t,«)) = f£ +1 (\Jm [Si]t n ^(i>, r , e )) c (Jf^ +1 (m [0]t n o^(r^ T , e )) , 

the image of f2jr(rjr T,e) under F^- +1 is a finite union of points belonging to M. Since, Vtjr{rjr T f ) is 
invariant, this a finite union of points, and thus a finite union of periodic orbits. 

The dynamics of neurons G T is driven by the periodic dynamics of firing neurons and it is easy to 
see that their trajectory is asymptotically periodic. Finally, since M. = U.y-lV,T,e the ui limit set of M is 
a finite union of periodic orbits. ■ 



A. 3 Average of a function. 

Since the dynamics is not uniquely ergodic (there are typically many periodic attractors), one has to 
be careful with the notion of average of a function (p. We have first to perform a time average for each 
attractor i, = fimT->oo Y^t=i ^(VW (t)), where is an initial condition in the attraction basin of 
attractor i. Then, we have to average over all attractors, with a weight corresponding to the Lebesgue 
measure /xW f its attraction basin. This gives: 

i=l 

where M is the number of attractors. 

References 

Abeles, M. (1982). Local cortical circuits: An electrophysiological study. Berlin: Springer. 
Abeles, M. (1991). Corticonics: Neural circuits of the cerebral cortex. New- York: Cambridge University 
Press. 

Abeles, M., Vaadia, E., Bergman, H., Prut, Y., Haalman, I., & Slovin, H. (1993). Dynamics of neuronal 

interactions in the frontal cortex of behaving monkeys. Concepts in Neuroscience, 4, 131-158. 
Amari, S. (1972, November). Characteristics of random nets of analog-like elements. IEEE Trans. Syst. 

Man and Cybernetics., SMC-2{5), 643-657. 
Ashwin, P., & Timme, M. (2005). Unstable attractors: existence and robustness in networks of oscillators 

with delayed pulse coupling. Nonlinearity , 18, 2035-2060. 
BenArous, C, & Guionnct, A. (1995). Large deviations for langevin spin glass dynamics. Probability 

Theory and Related Fields, 102, 455-509. 
Bertschinger, N., & Natschlagcr, T. (2004). Real-time computation at the edge of chaos in recurrent 

neural networks. Neural Computation, 16, 1413-1436. 
Blanchard, P., Cessac, B., & Krucgcr, T. (2000). What can one learn about self-organized criticality 

from dynamical system theory ? Journal of Statistical Physics, 98, 375-404. 



35 



Bressloff, P. C, & Coombcs, P. S. (2000a). A dynamical theory of spike train transitions in networks of 

integrate-and-fire oscillators. SI AM J Appl Math, 60, 820-841. 
Bressloff, P. C, & Coombes, P. S. (2000b). Dynamics of strongly coupled spiking neurons. Neural 

Computation, 12(1), 91-129. 
Brette, R., Rudolph, M., Carnevale, T., Hincs, M., Beeman, D., Bower, J. M., et al. (2007). Simula- 
tion of networks of spiking neurons: A review of tools and strategies. Journal of Computational 

Neuroscience, 23(3), 349-98. 
Broer, H., Efstathiou, K., & Subramanian, E. (2008). Robustness of unstable attractors in arbitrarily 

sized pulse-coupled systems with delay. Nonlinearity , 21(1), 13-49. 
Brunei, N., & Hakim, V. (1999). Fast global oscillations in networks of integrate-and-fire neurons with 

low firing rates. Neural Computation, 11, 1621-1671. 
Brunei, N., & Latham, O. (2003). Firing rate of noisy quadratic integrate-and-fire neurons. Neural 

Computation, 15, 2281-2306. 
Brunei, N., & Scrgi, S. (1998). Firing frequency of leaky integrate and fire neurons with synaptic current 

dynamics. J. Theor. Biol, 195(87-95). 
Cessac, B. (1994). Occurence of chaos and at line in random neural networks. Europhys. Lett., 26(8), 

577-582. 

Cessac, B. (1995). Increase in complexity in random neural networks. J. de Physique, 5, 409-432. 

Cessac, B. (2008). A discrete time neural network model with spiking neurons, rigorous results on the 
spontaneous dynamics. Journal of Mathematical Biology, 56(3), 311-345. 

Cessac, B., Blanchard, P., Krueger, T., & Meunier, J. (2004). Self-organized criticality and thermody- 
namic formalism. Journal of Statistical Physics, 115(516), 1283-1326. 

Cessac, B., Doyon, B., Quoy, M., & Samuclidcs, M. (1994). Mean-field equations, bifurcation map, and 
route to chaos in discrete time neural networks. Physica 74 D, 24-44. 

Cessac, B., & Samuelides, M. (2007). From neuron to neural networks dynamics. EPJ Special topics: 
Topics in Dynamical Neural Networks, 14-2(1), 7-88. 

Cessac, B., & Sepulchre, J. (2006). Transmitting a signal by amplitude modulation in a chaotic network. 
Chaos, ^(013104). 

Chow, C. C, & Kopell, N. (2000). Dynamics of spiking neurons with electrical coupling. Neural 

Computation, 12, 1643-1678. 
Coombcs, S. (1999a). Chaos in integrate-and-fire dynamical systems. In Proceedings of stochastic 

and chaotic dynamics in the lakes (Vol. 502, p. 88-93). American Institute of Physics Conference 

Proceedings. 

Coombes, S. (1999b). Liapunov exponents and mode-locked solutions for integrate-and-fire dynamical 

systems. Physics Letters A, 255, 49-57. 
Coombes, S., & Bressloff, P. C. (1999). Mode-locking and arnold tongues in integrate-and-fire neural 

oscillators. Physical Review E, 60, 2086-2096. 
Dayan, P., & Abbott, L. F. (2001). Theoretical neuroscience : Computational and mathematical modeling 

of neural systems. MIT Press. 
Derrida, B. (1987). Dynamical phase transition in non-symmetric spin glasses. J. Phys. A: Math. Gen., 

20, 721-725. 

Derrida, B., & Flyvbjerg, H. (1986). Multivalley structure in kauffman's model: analogy with spin 

glasses. J. Phys. A, 19, 1003-1008. 
Derrida, B., & Pomeau, Y. (1986). Random networks of automata: a simple annealed approximation. 

Europhys. Lett., 45-49. 



36 



Destexhe, A., & Contreras, D. (2006, October). Neuronal computations with stochastic network states. 
Science, 314, 85-90. 

Ernst, U., Pawelzik, K., & Geisel, T. (1995). Synchronization induced by temporal delays in pulse-coupled 

oscillators. Phys. Rev. Lett., 74{9), 1570-1573. 
FitzHugh, R. (1961). Impulses and physiological states in models of nerve membrane. Biophys. J., 1, 

445-466. 

Frcgnac, Y. (2004). From synaptic rumours to low-level perception: an intracellular view of visual cortical 

dynamics. Progress in Biochemistry and Biophysics, 31, 6-8. 
Frcidlin, M. I., & Wentzell, A. D. (1998). Random perturbations of dynamical systems. New York: 

Springer. 

Gambaudo, J., & Tresser, C. (1988). Le chaos, theorie et experiences. France: Collection CEA. 
Gerstner, W., & Kistler, W. (2002b). Spiking neuron models. Cambridge University Press. 
Gerstner, W., & Kistler, W. M. (2002a). Mathematical formulations of hebbian learning. Biol Cybern, 
87, 404-415. 

Golinelli, O., & Dcrrida, B. (1989). Barrier heights in the kauffman model. J. Physique, 50, 1587-1602. 
Golomb., D., Hertz, J., Panzeri, S., Treves, A., & Richmond, B. (1997). How well can we estimate the 

information carried in neuronal responses from limited samples? Neural Computation, 9, 649-665. 
Gong, P., & Leeuwen, C. van. (2007). Dynamically maintained spike timing sequences in networks of 

pulse-coupled oscillators with delays. Phys. Rev. Lett, 55(048104). 
Guckenheimer, J., & Holmes, P. (1983). Non linear oscillations, dynamical systems, and bifurcation of 

vector fields. Springer- Verlag. 
Hansel, D., & Mato, G. (2003). Asynchronous states and the emergence of synchrony in large networks 

of interacting excitatory and inhibitory neurons. Neural Computation, 15, 1-56. 
Hansel, D., Mato, G., Meunier, C, & Neltner, L. (1998). On numerical simulations of integrate-and-fire 

neural networks. Neural Computation, 10, 467-483. 
Hertz, J. (1997). Theoretical aspects of neural computation. In (p. 135-144). Wong K-Y M. King I. and 

Yeung D-Y (eds), Springer- Verlag. 
Hodgkin, A., & Huxley, A. (1952). A quantitative description of membrane current and its application 

to conduction and excitation in nerve. Journal of Physiology, 117, 500-544. 
Hunt, B. R., Saucr, T., & Yorkc, J. A. (1992). Prevalence: A translation-invariant 'almost every' on 

infinite-dimensional spaces. Bull. Am. Math. Soc, 27, 217-238. 
Jahnke, S., Mcmmcshcimcr, R.-M., & Timme, M. (2008). Stable irregular dynamics in complex neural 

networks. Phys. Rev. Lett., ^00(048102). 
Jolivet, R., Rauch, A., Lscher, H.-R., & Gerstner, W. (2006). "integrate-and-fire models with adaptation 

are good enough". MIT Press, Cambridge. 
Katok, A., & Hasselblatt, B. (1998). Introduction to the modern theory of dynamical systems. Kluwer. 
Kauffman, S. (1969). Metabolic stability and epigencsis in randomly constructed genetic nets. J. Theor. 

Biol, 22. 

Keener, J., Hoppensteadt, F., & Rinzel, J. (1981). Integrate-and-fire models of nerve membrane response 
to oscillatory input. SIAM Journal of Applied Mathematics, ^-Z(3), 503-517. 

Knight, B. (1972). Dynamics of encoding in a population of neurons. J. Gen. Physiol, 59, 734-766. 

Koch, C, & Segev, I. (Eds.). (1998). Methods in neuronal modeling: From ions to networks. The MIT 
Press. 

Kruglikov, A., & Rypdal, M. (2006a). Entropy via multiplicity. Discrete and continuous dynamical 
systems, 16(2). 



37 



Kruglikov, A., & Rypdal, M. (2006b). A piece- wise affine contracting map with positive entropy. Discrete 

and continuous dynamical systems, 16(2). 
Langton, C. (1990). Computation at the edge of chaos. Physica D., 42. 

Lapicque, L. (1907). Recherches quantitatives sur l'excitation electrique des nerfs traitee comme une 

polarisation. J. Physiol. Pathol. Gen., 9, 620-635. 
Latham, P., Roth, A., Hausser, M., & London, M. (2006). Requiem for the spike? Soc. Neurosc. Abstr., 

32. 

Maass, W., & Bishop, C. M. (Eds.). (2003). Pulsed neural networks. MIT Press. 

Mackay, R., & Tresser, C. (1986). Transition to topological chaos for circle maps. Physica D, 19(2), 
206-273. 

McCormick, D., Shu, Y., & Yu, Y. (2007, January). Neurophysiology: Hodgkin and huxley model - still 

standing? Nature, 445, E1-E2. 
Mcmmcshcimer, R.-M., & Timme, M. (2006). Designing the dynamics of spiking neural networks. Phys. 

Rev. Lett, 97(188101). 

Mezard, M., Parisi, G., & Virasoro, M. (1987). Spin-glass theory and beyond. World scientific Singapore. 
Mirollo, R. E., & Strogatz, S. H. (1990). Synchronization of pulse-coupled biological oscillators. SIAM 

J. Appl. Math., 50, 1645-1662. 
Muller, P. (1995). Calculation of lyapunov exponents for dynamic systems with discontinuities. Chaos, 

Solitons and Fractals, 5, 1671-1681. 
Nagumo, J., Arimoto, S., & Yoshizawa, S. (1962). An active pulse transmission line simulating nerve 

axon. Proc.IRE, 50, 2061-2070. 
Naundorf, B., Wolf, F., & Volgushev, M. (2006, April). Unique features of action potential initiation in 

cortical neurons. Nature, 440, 1060-1063. 
Naundorf, B., Wolf, F., & Volgushev, M. (2007, January). Neurophysiology: Hodgkin and huxley model 

- still standing? (reply). Nature, 445, E2-E3. 
Packard, N. (1988). Dynamic patterns in complex systems. In (p. 293-301). Kelso J.A.S, Mandcll A.J. 

and Shlesinger M.F. Editors, World Scientific. 
Panzeri, S., & Treves, A. (1996). Analytical estimates of limited sampling biases in different information 

measures. Networks, 7, 87-107. 
Rieke, F., Warland, D., Steveninck, R. de Ruyter von, & Bialck, W. (1996). Spikes, exploring the neural 

code. The M.I.T. Press. 

Rudolph, M., & Destexhe, A. (2006). Analytical integrate and fire neuron models with conductance-based 
dynamics for event driven simulation strategies. Neural Computation, 18, 2146-2210. 

Samuelides, M., & Cessac, B. (2007). Random recurrent neural networks. EPJ Special Topics "Topics in 
Dynamical Neural Networks : From Large Scale Neural Networks to Motor Control and Vision", 
142(1), 89-122. 

Senn, W., & Urbanczik, R. (2001). Similar non-leaky integrate-and-fire neurons with instantaneous 

couplings always synchronize. SIAM J. Appl. Math., 61(A), 1143-1155. 
Siri, B., Berry, H., Cessac, B., Delord, B., & Quoy, M. (2007). Effects of hebbian learning on the dynamics 

and structure of random networks with inhibitory and excitatory neurons. Journal of Physiology 

(Paris), 101 ((1-3)), 138-150. (c-print: arXiv:0706.2602) 
Siri, B., Berry, H., Cessac, B., Delord, B., & Quoy, M. (2008). A mathematical analysis of the effects 

of hebbian learning rules on the dynamics and structure of discrete-time random recurrent neural 

networks. Neural Computation, (to appear) 
Sompolinsky, H., Crisanti, A., & Sommcrs, H. (1988). Chaos in random neural networks. Phys. Rev. 

Lett, 61, 259-262. 



38 



Soula, H. (2005). Dynamique et plasticite dans les reseaux de neurones a impulsions. Unpublished 

doctoral dissertation, INSA Lyon. 
Soula, H., Beslon, G., & Mazet, O. (2006). Spontaneous dynamics of asymmetric random recurrent 

spiking neural networks. Neural Computation, 18(1). 
Timme, M., Wolf, F., & Geisel, T. (2002). Coexistence of regular and irregular dynamics in complex 

networks of pulse-coupled oscillators. Phys. Rev. Lett., 89, 258701. 
Touboul, J., & Faugeras, O. (2007). The spikes trains probability distributions: a stochastic calculus 

approach. Journal of Physiology, Paris. (To appear.) 
VanVreeswijk, C. (2004). What is the neural code? 23 Problems in System neuroscience. van Hemmcn, 

J.L. and Sejnowski, T.Jr. (cds), Oxford University Press. 
VanVreeswijk, C, & Hansel, D. (1997). Rhythmic bursting in networks of adaptive spiking neurons. 

Computational Neuroscience 97 Abstracts. 
VanVreeswijk, C, & Sompolinsky, H. (1998). Chaotic balanced state in a model of cortical circuits. 

Neural Computation, 10, 1321-1372. 



39 



